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Abstract 

We present a procedure that in many cases enables the Monte Carlo sampling of 
states of a large system from the sampling of states of a smaller system. We illustrate 
this procedure, which we call the sewing algorithm, for sampling states from the 
transfer matrix of the two-dimensional Ising model. 
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1 Introduction 



From its onset, the Monte Carlo method was recognized as an effective ap- 
proach for estimating the solutions to linear systems of equations, inverting a 
matrix, and finding eigenvalues of a matrix when these linear algebra problems 
involved a very large matrix PQ . The core operation in estimating the solutions 
of such problems is performing a sequence of matrix-vector multiplications, 

^ ^ -A-ijAjk ' * ' A mn x n . (1) 

j,k....,n 



Instead of executing it completely, the Monte Carlo approach samples terms in 
the summations, for example, A ili2 A i2i3 ■ ■ ■ A in _ lin Xi n , by generating a Markov 
chain (zi, . . . , i n ) denned by a matrix Tij for transition from state j to i and 
relating the random walks generated by to the sequence matrix- vector mul- 
tiplications by having the walker accumulate a "weight" Wi^Wi^ ■ ■ ■Wi n _ x i n 
as it moves along the chain. From the samples, expectation values of the result 
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are easily constructed. There is considerable leeway in choosing the Tij. The 
minimum conditions are 

Aij = TijWij (2) 

where is greater than zero if Ay is nonzero, equals zero if does, and 
satisfies X^^V/ = 1 for all j. The core Monte Carlo sampling is for a given j 
selecting an i with the conditional probability T^. 

Sampling the conditional probability is straightforward until the matrix be- 
comes too large to store or too expensive to continually regenerate. It is these 
situations that this paper addresses by presenting a new strategy, the sewing 
algorithm, which should be applicable to a variety of problems. The sewing al- 
gorithm accomplishes the sampling of states i for a large system from samples 
of states from transition matrices of small systems. It was developed in the 
context of a power method algorithm to obtain multiple eigenvalues of very 
large matrices |2,3,4.5j. In this eigenvalue work the benchmark application 
was obtaining the two largest eigenvalues of the transfer matrix of the two- 
dimensional Ising model [6]. With the Sewing Algorithm we easily obtained 
these two eigenvalues for a 60 x 60 lattice on a single processor. Past work 
required parallel computing. 60 x 60 is where we decided to stop and likely not 
the limit of the algorithm. Reference [5] contains details about the application 
and references to the eigenvalue algorithms. 

We will present the sewing algorithm in the context of the transfer matrix for 
the two-dimensional Ising model. While some of the details will be relevant 
only to this and related problems, the basic strategy is more general. In the 
next section, Section 2, we will define the transfer matrix [TJ. In Section 3, we 
will discuss the basics for implementing the Monte Carlo to sample a sequence 
of matrix-vector multiplications when the system is still small enough so that 
there is enough computer memory to store T^. Then, in Section 4, we will 
discuss how to sample T^, when memory is inadequate for the system size of 
interest but is adequate for smaller system sizes. Finally, in the last section, 
we will make some comments about the generality of the method. 



2 Transfer Matrix of Two-Dimensional Ising Model 

We will consider an m x m Ising model defined with periodic boundary condi- 
tions in one direction and open boundary conditions in the other. The model's 
energy is [6] 

m— 1 m mm 

E M = -J J2 J2 - J J2J2 IMj&j+i ( 3 ) 

i=l j=l i=l j=l 



2 



Here, are the coordinates of a lattice site. The Ising spin variable /i^j on 
each site has the value of ±1, the exchange constant J > 0, and /ij jm +i = 
The symbol 



(4) 



denotes a column configuration of Ising spins and there are 2 m possible config- 
urations for each column. Typically, a configuration is mapped onto an integer 
in the range to 2 m — 1 with a +1 Ising spin mapped to a 1 bit in the in- 
teger and a —1 spin to a bit. The transfer matrix A(a, o') follows from a 
re-expression of the partition function [6J 



Z (m, m)=Y^ exp [-vE ({//})] 



ex P v wl {^i + 5 2 (Oi, t7 i+ i)} 
o"i,...,<T m Vi=i 

: 51 -4(o-i, cr 2 )A(o- 2 ,o-3) • • • A(cr m _i,a m )A(cr m ,ai; 



(5) 



where 

^1 ( a j) = Vij^i+lj 



m—l 



i=l 



is the interaction energy of the j th column and 

m 

Mi j'/^i J+l 



S 2 (aj,a j+1 ) 



i=i 



(6) 



(7) 



is the interaction energy between the j th and (j + l) th columns, v = J/ksT 
kg is Boltzmann's constant, and T is the temperature. A(a, a') is a 2 m x 2" 
matrix whose elements are 



m—l 



A (o, a') = exp v ^ /ifc/^+i exp v ^ 



fc=i 



fc=i 



Because of the one open boundary, the matrix is asymmetric. We note that 
all the elements of A(a, a') are greater than zero. The mapping of the o to 
integers maps A (o, a') to A^. 
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3 Core Monte Carlo 



We assume that the elements of some MxM matrix A are easily generated on- 
t he-fry. For simplicity, but without loss of generality, we will also assume that 
all its elements are positive. Next, we imagine we have N random walkers 
distributed over the M states defining A, and we will interpret Ay as the 
weight of particles arriving in state i per unit weight of a walker in state j and 
will regard the action of A on x as causing a walker to jump from some j to 
some i, carrying its current weight Xj, modified by Aij, to state i. Repeated 
walks by a walker and the sum over all walkers estimates ([T]) . 

The jumps are executed probabilistically. To do this, instead of (J2J), we define 
the total weight leaving state j as 

^i = E4* (9) 

i 

and the transition probability from j to % as 

T tj = A lJ /W j (10) 



The number Wj is called the state weight multiplier. 

If M is sufficiently small, then a Monte Carlo procedure is easily constructed. 
We sample a state i from T^- and multiply the transferred weight by the ratio 
of the true probability (1.0) to the sampled probability (T^); that is, if state 
i is sampled, the weight arriving in state i from state j is multiplied by 



The standard way [8] to sample % from T is to first construct the cumulative 
distribution, 

i 

C ij = "l2 T kj (12) 
k=l 

and then draw a random number £ from the unit distribution. The state i 
equals the value of k that satisfies 

Cfc-i,j < £ < Ckj (13) 



where Cqj = 0. 
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4 Sewing Algorithm 



For M small enough so can be stored in memory, sampling from the 
cumulative probability CV, works well. If the number of states gets too large, 
then Cij cannot be sampled directly In this case, instead of, for example, 
just randomly picking from any state j any state % with probability 1/M, 
we developed a new procedure which we call the sewing algorithm. In it, we 
assume, for example, that we can write any state iasa direct product of the 
states in a smaller basis. If it is the direct product of two states, that is, 

\i) = |i 2 > \h) (14) 



then instead of transferring weight 

w,- = £4y (15) 

k 



from state j to state % with probability 

T ii =A ij / , £Ak j = A ij /W jl (16) 
k 



we will use the a^- that would apply to the smaller set of states and then make 
an appropriate weight correction. 

For each smaller set of states, we have 



t ij = a ij /w j , (17) 



and for \j) = \j 2 ) \ji) we sample and \i 2 ) from 

(19) 

Now, we define to be the weight correction necessary to preserve the 
expected weight transfer from state j to state i. It satisfies 

Aij = Xijti 2 j 2 ti 1 j 1 (20) 
Thus 

A ii = X ij ai2j2ailjl (21) 
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that is 

A- ■ 

X i: j = w h w j2 %1 — ( 22 ) 

a hj 1 a i 2 j2 



This sewing method generalizes easily. For k sets of states, ff2TT) and ff22l 
become 

k 

Aij = Y[ U n j„ (23) 

n=l 

and 

k in ■ 

*ij = Mi II (24) 



4-1 Illustration for the Transfer Matrix 

For the Ising model 

A ij = esxp[v{S 1 {i)+S 2 {i,j))] (25) 

where S± and S 2 were defined in and ([7j). To compute the weight correction 
factor X^, we need to study how the sums Si(i) and S 2 (i,j) differ between 
evaluating them with the bits taken together and taken separately. 

We will establish this difference in the context of the following assumed binary 
representation of two states 

\i) = 110101000110011011110) = |1010100011)|0011011110) (26) 
\j) = 110101100111011011110) = |1010110011)|1011011110) (27) 

The sum Si(i), computed on i alone, adds 1 every time adjacent bits match 
and subtracts 1 every time they mismatch. Because of the periodic boundary 
condition, the two ends are adjacent. For ([26"]) . there are 12 matches and 8 
mismatches so S\{i) = —4. The sum S 2 (i,j) counts matches and mismatches 
between the bits of i and j. For (f26|) and (l27j) . there are 18 matches and 2 
mismatches so S 2 (i,j) = 16. The transfer matrix element from j to % is then 

= expKS^) + S 2 {i,j))\ = exp(12z/) (28) 

In sewing the bits together, our intent was to proceed recursively so we asumed 
periodic boundary conditions to hold for the various bit segments; that is, 
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was the transfer matrix of a smaller-sized Ising model. First, we will consider 
just the lower 10 bits of i. 



There are 6 matches and 4 mismatches so its first sum si(zi) = 2. For the 
upper 10 bits of i, 



there are 4 matches and 6 mismatches so its first first sum 51(22) = —2 

We note that the terms making up Si(i) are almost the same as the terms 
making up sj.(zi) and sifa)- The differences occur at the ends of the sets 
of bits. In the present example the difference between taking bit segments 
together rather than separately are two mismatches occurring at the right 
ends of each segment. Thus Si(i) = (si(ii) — 2) + (si(i 2 ) — 2) = —4. We now 
write Si(i) = Si(i\) + si(i 2 ) + D(i) where D(i) is defined to be the "energy" 
difference between calculating the sets of bits together and calculating the sets 
of bits separately. 

For S 2 (i,j) the number of matches and mismatches between i and j is just 
the sum of the matches and mismatches on the sewn pieces, that is, S 2 (i, j) = 
s 2 («i,ji) + s 2 (i 2 , jz). 

We can now write 



ii) = 10011011110) 



(29) 



i 2 ) = 11010100011) 



(30) 



X i: j = w j2 w jl 



exp[u(Si{i) + S 2 (i,j))} 



(31) 



exp[i/(si(z 2 ) + s 2 (i 2 ,j 2 ))] exp[i/(si(z'i) + s 2 (h, ji))} 



We showed that 



Sz{i,j) = s 2 (i 2 ,j 2 ) + s 2 (h, ji) 



(32) 



and 



Si(i) = si(i 2 ) + si(zi) + D(i) 



(33) 



Substituting (|32j) and (!33|) into (I3TT) yields 



Xij = Wj^Wj^ exp [vD{i)] 



(34) 
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The weight correction for sewing k systems together becomes 

k 

Xij = exp(uDi) Yl Wj n (35) 

n=l 

where vDi is the energy difference per ksT between calculating with the bits 
together and the bits separately. 



5 Concluding Remarks 

The sewing algorithm is very simple and precise in concept, and its application 
should extend beyond the specific case illustrated. Our presentation had two 
classes of assumptions, one for the matrices and the other for the states. 

We assumed all our matrix elements were positive. If some are zero, then 
the corresponding T^- must be zero. If some are negative, the basic strategy 
prevails when absolute values are taken in the appropriate places. Mixed signed 
matrix elements necessitates mixed signs of the weights of the walkers, and 
additional Monte Carlo procedures might be necessary to promote proper 
cancellations of the signs. This need was the case for our work in computing 
multiple eigenvalues of large matrix via the Monte Carlo method [5] . 

We also assumed the state of the larger system was expressible as a direct 
product of the states of a smaller systems. Clearly, this specific requirement 
was possible and convenient for the Ising problem and is likely so for other 
problems involving interacting spins, electrons, and hard-core bosons. More 
generally, the larger system simply needs to be domain-like decomposable 
with interfacial corrections easily computed. The details for doing this will 
depend of the specific problem. Many, but not all, problems will be efficiently 
amendable to the sewing algorithm. 
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